############################## 森林图 #############################
##### 科研工具8 森林图
library(missRanger) #用于缺失值处理
library(haven) #用于读取xpt
library(gtsummary) # 统计
library(forestplot)
# 测试年龄与体重的关系 协变量 性别 收入 学历
# 年龄 RIDAGEYR
# 性别 RIAGENDR
# 教育水平 DMDEDUC2
# 收入 INDFMPIR
getwd()
setwd('父级路径')
# ismy <- file.exists("forestest.csv")
# if(ismy){
#  
#   setwd("F:/rproject/pmwosplot")
# }else{
#   setwd("G:/BaiduNetdiskDownload")
#   demo.d <- read_xpt("NHANES/2005-2006/Demographics/demo_d.xpt")#参见上述设置默认路径，文件名称后缀不同"d,e"
#   demo.e <- read_xpt("NHANES/2007-2008/Demographics/demo_e.xpt")#参见上述设置默认路径，文件名称后缀不同"d,e"
#   bmx.d <- read_xpt("NHANES/2005-2006/Examination/bmx_d.xpt")#参见上述设置默认路径，文件名称后缀不同"d,e"
#   bmx.e <- read_xpt("NHANES/2007-2008/Examination/bmx_e.xpt")#参见上述设置默认路径，文件名称后缀不同"d,e"
#   #合并两个cycle 数据
#   demo_data_all <- dplyr::bind_rows(list(demo.d,demo.e))
#   bmx_data_all <- dplyr::bind_rows(list(bmx.d,bmx.e))
#   demo_data <- demo_data_all[,c('SEQN', 'RIDAGEYR', 'RIAGENDR', 'DMDEDUC2', 'INDFMPIR'  )]
#   bmx_data <- bmx_data_all[,c('SEQN', 'BMXBMI'  )]
#   output <- plyr::join_all(list(demo_data, bmx_data), by='SEQN', type='full')
#   # 清除掉空值的数据
#   output <- output[which(!is.na(output$BMXBMI)),]
#   setwd('F:/Rproject/r-language/paper1Somke/')
#   write.csv(output,'forestest.csv')
# }

# agerange317dropna <- read.csv('source.csv')

output <-read.csv('forestest.csv')


colnames(output)


# 默认线性回归模型
m1 <- glm(BMXBMI ~ RIDAGEYR + RIAGENDR+DMDEDUC2+INDFMPIR, data = output)
#下面是逻辑回归模型 只是加了一个参数
#m12 <- glm(BMXBMI ~ RIDAGEYR + RIAGENDR+DMDEDUC2+INDFMPIR, data = output,family = 'binomial')

reg <- tbl_regression(m1,exponentiate =T )
#reg2 <- tbl_regression(m2,exponentiate =T )

try({
  reg%>%as_flex_table()%>%flextable::save_as_image(path = 'glmpic.tiff')
})


# 服务器报错 后续找原因


try({
  gt::gtsave(as_gt(reg), 'glm.png')
})


# gt::gtsave(as_gt(reg), file = file.path(getwd(), "glm.png"))

# Cairo::CairoTIFF(file="foresttest.tiff", width=8, height=8,units="in",dpi=150)
# plot(reg)
# dev.off()

fit.result<-summary(m1)
df1<-fit.result$coefficients
df2<-confint(m1)
df3<-cbind(df1,df2)
df4<-data.frame(df3[-1,c(1,4,5,6)])
df4$Var<-rownames(df4)
colnames(df4)<-c("OR","Pvalue","OR_1","OR_2","Var")
df5<-df4[,c(5,1,2,3,4)]
df5$OR_mean<-df5$OR
df5$OR<-paste0(round(df5$OR,2),
               "(",
               round(df5$OR_1,2),
               "~",
               round(df5$OR_2,2),
               ")")
df5$Pvalue<-round(df5$Pvalue,3)


Cairo::CairoTIFF(file="BMIforest.tiff", width=8, height=8,units="in",dpi=150)
forestplot(labeltext=as.matrix(df5[,1:3]),
           mean=df5$OR_mean,
           lower=df5$OR_1,
           upper=df5$OR_2,
           zero=0.25,
           boxsize=0.2,
           graph.pos=2)
dev.off()

